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We introduce a Monte-Carlo algorithm for the simulation of charged particles moving in the continuum. 
Electrostatic interactions are not instantaneous as in conventional approaches, but are mediated by a constrained, 
diffusing electric field on an interpolating lattice. We discuss the theoretical justifications of the algorithm and 
show that it efficiently equilibrates model polyelectrolytes and polar fluids. In order to reduce lattice artifacts 
that arise from the interpolation of charges to the grid we implement a local, dynamic subtraction algorithm. 
This dynamic scheme is completely general and can also be used with other Coulomb codes, such as multigrid 
based methods. 
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I. INTRODUCTION 

Simulating charged condensed matter systems is demand- 
ing due to the long range nature of the Coulomb interaction'-^. 
The direct evaluation of the Coulomb sum for N parti- 
cles, Uc = ^i<j ejCj"/47reorij, requires computation of 
the separations between all pairs of particles, which im- 
plies ©(iV^) operations are needed per sweep or time step. 
Conventional fast algorithms including the classic Ewald 
sum.^, particle mesh Ewald approaches'*"!, and fast multipole 
algorithms^ suffer from either poor scaling with system size, 
high coding complexity or inefficiency in multiprocessor en- 
vironments. Although multigrid methods^ appear to be a 
promising candidate for fast parallelized molecular dynam- 
ics simulations, these techniques are less efficient for Monte- 
Carlo simulations, where the total energy of the system must 
be recalculated after every particle move. 

It is surprising that the Coulomb interaction poses such 
tremendous difficulty; after all the underlying Maxwell equa- 
tions are local. The above methods, however, do not solve 
Maxwell's equations, but rather search for the electrostatic po- 
tential (j)p{r), from which the electric field E(r) is deduced. 
We shall show in this article that this full relaxation of the 
electric field to the solution of the Poisson equation is not 
necessary in the context of statistical mechanics, where one is 
interested in configurational averages. This allows us to con- 
struct a nonstandard approach that is local and substantially 
more efficient than methods involving potential theory. 

Recently, a lattice Monte-Carlo algorithm was introduced 
in which the time needed for a single sweep scales as only 
0{N). This algorithm employs the electric field, E, rather 
than the scalar potential as the true degree of freedom.**'*. 
For a charge density p{r) = J^i eiS{r — r^). Gauss' law, 
eodivE(r) — p{r) = 0, is imposed as an exact dynamical 
constraint on the field configurations. The algorithm correctly 
reproduced the static properties of a charged lattice gas and 
equilibrated the system with a 0(1) overhead when compared 
to the time needed to simulate a neutral system.'". This high 
efficiency results from purely local updates of particles and 
fields using the standard Metropolis Monte-Carlo rule; there 
is no global calculation of the potential. 

In the present article, we generalize this algorithm to off- 
lattice Monte-Carlo, so that the particles move in the contin- 



uum. Such a generalization is necessary for any realistic mod- 
eling of charged systems, such as polyelectrolytes, biopoly- 
mers and polar solvents. We begin with a complete descrip- 
tion of the algorithm in Section HI] and continue with tests of 
the static and dynamic behavior in Section In particular, 
we demonstrate numerically that the slowest relaxation time 
in the charged system is dominated by diffusion of the den- 
sity degrees of freedom, in a manner which is very similar to 
a neutral system with short range interactions. The following 
Section IIVI discusses the elimination of lattice artifacts that 
arise from the interpolation of charges to a grid via dynamic 
subtraction. Although demonstrated in the context of our aux- 
iliary field Monte-Carlo algorithm, this dynamic subtraction 
can also be used with more conventional Coulomb algorithms 
such as multigrid. 



n. COULOMB INTERACTIONS FROM AUXILIARY 
FIELDS 

A. Constrained energy functional 

The starting point for our derivation of a local Coulomb 
algorithm is a formulation of electrostatics in terms of a con- 
strained energy functional based on a vector field E(r)ii, 



^fEl = 



eoE(r)2 



-0(r)(eodivE(r)-p(r)) 



d^r. (1) 



In this functional Gauss' law is imposed with the Lagrange 
multiplier (j> and constrains the longitudinal component of the 
electric field; 0(r) can be identified with the scalar potential 
(f>p. Minimization of the functional with respect to E leads 
immediately to the familiar equations of electrostatics. 
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All conventional Coulomb algorithms are essentially con- 
cerned with finding solutions to Eqs. ( I2I3> . After inserting 
Eqs. ( I2I3> into Eq. Q, we obtain the electrostatic energy 
Uc — Y j(S4>p)'^d^''' as a minimum of the constrained func- 
tional. Note that a scalar functional based only on the potential 



2 



(j) cannot be used, it leads to the wrong sign in the Coulomb 
interaction (see Ref. 10 for a discussion of this point). 

The crucial point of our algorithm is now not to minimize 
the constrained energy functional Eq. Q, but rather to main- 
tain the electric field at finite temperature. For T > 0, we have 

E = -V(/)p + Et^, (4) 

where E^^ = curlQ is an arbitrary rotational vector 
field. The total field E still satisfies Gauss' law, but 
the energy U is no longer equal to Uc, but rather U = 
^ J {(V^ip)^ + Ej^} d^r. We now show that an integration 
over these auxiliary transverse (rotational) degrees of freedom 
of the electric field nevertheless lead to electrostatic interac- 
tions. 

To this end, let us consider the partial partition function for 
charged particles and a fluctuating electric field that is con- 
strained by Gauss' law, 

Z({rJ) - /pE n'^(divE-p({r,})/eo)e- ^'''"''«^'/''-«^. 

r 

(5) 

Note that the integration is only performed over field configu- 
rations and not particle positions. We evaluate the integral by 
changing variables, Et^ = E + V(/)p, 

Z({rJ) = /"pEt,[|j(divEt,)e-/^''^''°(^"-^'^-)'/2fe^^ 

•' r 

^ p-Jd^r'eo(V0p)V2feBT 

X / VEtr n '^(div Etr) e- ^ ^o^J^k^T 

•' r 

= Zcouio7nb{{ri}) X const. (6) 

The cross term in the energy, J Et^ V^pd^r, vanishes for peri- 
odic or constant potential boundary conditions, as can be seen 
by integration by parts. Eq. (|6} shows that an integration over 
the transverse modes merely multiplies the coulombic parti- 
tion function Zcouiomb by a constant. Since this multiplica- 
tion does not change the relative statistical weight of configu- 
rations, it can be ignored. 

This result allows us to maintain the total electric field at 
finite temperature rather than quenching it to zero as when 
solving Poisson's equation. As long as the constraint of 
Coulomb's law is strictly maintained, an integration over all 
possible transverse field components will generate an effec- 
tive Coulomb interaction between the particles. As illustrated 
below, this integration involves only local operations and re- 
quires no global optimization. Although this may sound un- 
usual, it is actually a very natural way of calculating the inter- 
actions: we shall see below that the resulting dynamics has a 
local conservation law exactly as occurs in Maxwell's equa- 
tions. 



B. Algorithm 

Our Monte-Carlo algorithm evaluates the partition function 
for a Coulomb system using the Metropolis criterion together 



with the energy functional Eq. Q. The algorithm consists of 
two Monte-Carlo moves that 

• integrate over the particle positions by displacing the 
charged particles, but maintain a field configuration that 
always satisfies Gauss' law; 

• integrate over the transverse components of the electric 
field Eitr to evaluate the configurational integral Eq. (|6|l. 

We shall now describe the implementation of the algorithm 
for the off-lattice case. A similar description of the algorithm 
for lattice simulations has already been given in Ref. 10. The 
reader who is only interested in the results of the algorithm 
may jump directly to Sectionlllll 

1. Discretization 

Particles with charge e.i — ±1 move continuously in a cu- 
bic, periodic simulation cell of volume L^. The electric field 
E is discretized onto a cubic grid of mesh size a, and the 
Cartesian components of E are associated with the 3L^/a^ 
links of the lattice. A natural unit of temperature or energy 
is T* = / A-Ke^ksa. In these units, the Bjerrum length 
is l\j = aT* /T. A typical value in water at room tempera- 
ture is 4 ~ 7 A for monovalent ions. Charges are interpo- 
lated onto the /a^ nodes of the lattice. Various choices 
are possible for the interpolation scheme: we desire a rea- 
sonably smooth interpolation and choose a 3rd order scheme 
using B-splines^ '^. In each Cartesian direction, the charge is 
distributed onto the 3 mesh points closest to the particle. Ex- 
phcitly, the one-dimensional weights are 

W^{[x]) = (N~0.5)V2 
WoM) = 0.75- 

W+{[x]) = ([.t] + 0.5)V2, (7) 

where [x] denotes the distance between the particle and the 
nearest lattice point in units of the lattice spacing. The three- 
dimensional weights are constructed by multiplying the three 
one-dimensional weights. Our method is not tied to a particu- 
lar assignment scheme; higher or lower order schemes can be 
employed as long as they conserve charge exactly. 

2. Particle motion and Hamiltonian paths 

The simulation begins with a field configuration that satis- 
fies the discretized, integral version of Gauss' law, 

i 

In this expression Ei,j = —Ej^i denotes a field component 
leaving site i towards site j, and the sum is performed over 
all links connecting to site i; it corresponds to the total flux 
leaving site i. The success of the algorithm depends on main- 
taining the constraint at all times; the field configuration is 
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J= Aq^/a^ 
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FIG. 1: Illustration of charge interpolation and the current update 
using a Hamiltonian path in two dimensions. A charge (open circle) 
interpolates onto the 3^ = 9 lattice sites (solid). The partial charges 
on these sites change by an amount Aqi when the particle moves. 
The current J flows along the Hamiltonian path and terminates at the 
last site, since J2i ~ 0- 



updated each time a particle moves. For a lattice gas, such 
an update is easy to perform.*^; when a particle moves from 
site i to site j, one changes the field on the traversed link to 
Ei,j Eij ^ e/ {e^a?) (the sign depends on the direction in 
which the particle moves). The resulting field configuration 
once again satisfies Eq. (|8|l. The situation is more compli- 
cated in the off-lattice case; a single particle move changes the 
charges on 3'^ = 27 nodes, more if the particle changes cells. 
One can choose among many update rules to realign the field 
configuration at all modified nodes (the system of equations 
is in general under-specified). This task can be accomplished 
with minimal effort using a path that traverses the modified 
region in such a way that it visits each site once. An example 
of such a path, which is often called a "Hamiltonian path", 
is shown for two dimensions in Fig. ^ the generalization to 
three dimensions is obvious. Only links along the path need 
be updated. 

The field Ei,j on a link that is part of the path is updated 
according to Eij Ei^j ± Aq^/ (eoa^), where Aqs de- 
notes the change of the fractional charge on a node due to the 
particle motion. The primed summation extends over all sites 
that were previously visited by the path to reach site i, and 
the sign depends on the direction of the path. This update is 
efficiently performed during a single "walk" along the path. 
In a continuum description, the current that flows due to the 
motion of the particles produces a change in the electric field. 




FIG. 2: Integration over the transverse degrees of freedom. The 
four link fields forming a plaquette in the xy-plane are modified 
by ±A in such a way that the total field changes by a pure ro- 
tation (curlE)z = 4A. A is chosen from a uniform distribu- 
tion. We associate an angular variable Oz with the plaquette so that 
5E = curl<5©. 



Note that though the path changes direction several times dur- 
ing the update, the resulting currents are essentially parallel. 
We always construct the Hamiltonian path in a manner which 
is symmetric between the new and old interpolations sites. 
Because of this our procedure is time reversible and the al- 
gorithm respects detailed balance. To simplify the implemen- 
tation of the algorithm we always displace particles parallel 
to the links of the mesh with a step size uniformly distributed 
between and a. 



3. Integration over the transverse field 

Field updates as described in the previous paragraph change 
not only the longitudinal component of the electric field, but 
simultaneously affect the transverse field. The particle motion 
thus already performs a partial integration over the transverse 
modes. A complete evaluation of the configurational integral 
Eq. (|6j often requires a second type of Monte-Carlo move that 
changes just the transverse degrees of freedom of E without 
moving a particle nor violating the constraint Eq. (jS). This is 
done by grouping the links that form a faces of the grid into 
plaquettes. We randomly choose a plaquette and incre- 
ment the field on two of the links by an increment A while 
decreasing it on the other two; we only modify the total field 
configuration by a pure rotation (see Fig.|2j. We describe the 
resulting dynamics by associating an angular variable 9k with 
a plaquette in the {i, j}-plane. As one can see by inspection 
of Fig.|2] changes in the electric field are related to the curl of 
the changes of the plaquette variable 9k- Combining this with 
Eq. the evolution of the electric field obeys the equation 



'dt 



-J/eo - 



(10) 



9E 

'dt 



-J/eo. 



Changes in E also couple back to the circulatory degrees of 

(9) 

freedom, one can show^ that curl E acts like a torque on 9, i.e. 
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^— = -eocurlE + 77, 



(11) 



where we have introduced a relaxation time scale ^ which is 
determined by the Monte-Carlo moves; 77 is a Brownian noise. 

The dynamics of Eqs. JlOll 1> is closely related to the 
Maxwell equations". After identifying 89/ dt with the mag- 
netic field B//io> we recover Ampere's law in Eq. ilO\ . If we 
perform the same replacement in Eq. il 1> . we recognize its 
closeness to another Maxwell equation, Faraday's law. How- 
ever, the usual ffR/dt has been replaced by just B; a time 
derivative has been lost. By combining the two equations one 
findsS that the electric field obeys a diffusion equation 



5E 

— = (eoV^E - Vp)/e - J/eo + curl^. 
ot 



(12) 



Since in Monte-Carlo we are interested in distributions rather 
than dynamics we are free to modify the underlying dynam- 
ics. Gauss' law is imposed as an initial condition and locally 
conserved both in Maxwell's equations and in our algorithm, 
as can be seen by taking the divergence of Eq. M2\ . The con- 
servation of Gauss' law in our algorithm and in Maxwell's 
equation is the origin of Coulomb's law. 



C. Boundary Conditions 

Traditional treatments of the Coulomb interaction in pe- 
riodic systems lead to a residual dependence on the dielec- 
tric constant e' of a surrounding medium.'^ ''*. For a standard 
Ewald summation, the short-range expansion of the pair po- 
tential between two point charges of opposite sign separated 
by r iai^ 



V{r) 



1 



(13) 



The present algorithm imposes a different choice that we de- 
note "Maxwell" boundary conditions, since they are the same 
as those generated by the full Maxwell equations. We show in 
the appendix that 



V{r) 



47reQ 



0(r4) 



for 
for 



(14) 



There is again a quadratic correction to the bare 1/r inter- 
action, which, surprisingly, depends on the temperature. A 
crossover occurs when the Bjerrum length 4 reaches the sys- 
tem size L. It is possible to remove the temperature depen- 
dence of the potential by introducing a third type of Monte- 
Carlo move that couples to the mean or q = compo- 
nent of the electric field E = -p- J <Pr'E,{r). The potential 
then coincides with the conventional Ewald sum with "tinfoil" 
(e' = co) boundary conditions. 
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FIG. 3: Pair potentials obtained from the distribution P(r) at dis- 
tance r of a pair of charges for two system sizes (a) L = 5 a and (b) 
L = 10 a and two different temperatures T; = O.OGST* (□, ★) and 
Th = IT* (A, 0). Solid lines Eq. (14}, V{r) = const - l/47rr- - 
h{T)r'^, where the limiting values for the quadratic correction are 
h{Ti) = -l/3eoi^ and h{Th) = l/6eoL^ for Maxwell boundary 
conditions (□, A) (dashed lines show the curves in the opposite limit 
to demonstrate sensitivity). Addition of a third move leads to tinfoil 
boundary conditions (*, 0), b = — l/3eoi^ for all T. All curves 
end at r = 2^^^' a, where the particles interact with the truncated 
Lennard- Jones potential. 



ni. RESULTS 
A. Static quantities 

1. Pair potential 

We begin the analysis of the algorithm by verifying that the 
effective interaction between particles is in agreement with 
Eq. (I14> . We evaluated the pair potential by placing two par- 
ticles of opposite charge in the simulation cell and measuring 
the distribution P{r) of separations between the particles for 
different temperatures and system sizes. The negative loga- 
rithm of P{r), which is proportional to the pair potential, is 
shown in Fig. |5] together with the curve Eq. J14> . To include 
self-avoidance, a purely repulsive truncated Lennard-Jones 
potential of range a was added to the interaction. Overall, 
we find excellent agreement between Eq. ( I14> and the numer- 
ical curves. For Maxwell boundary conditions, the quadratic 
correction depends on temperature, and the two temperatures 
shown in Fig.|3lcorrespond to the limiting values of Eq. ( I14> . 
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FIG. 4: Static charge-charge structure factor for a simple Coulomb 
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0.1 and p = 0.2 a'^ The slopes at large values of l/q^ in- i ^. exp (iq • Vi{t)), where the weights are replaced 



B. Efficiency 

Having shown that the algorithm behaves correctly, we now 
turn to a study of the dynamics. As already discussed in^, the 
Monte-Carlo update rules lead to a diffusive dynamics of the 
longitudinal and transverse fields, and we wish to explore how 
quickly our off-lattice algorithm equilibrates various modes of 
the system. 



I. Dynamics of free charges 

We investigate the dynamics of the simple Coulomb gas 
shown in Fig. 0] by studying the time decay of correla- 
tions in Fourier space. We measure the dynamical struc- 
ture factors S'(q, t) = (s(q, t)s(— q, 0)) with s(q, i) = 



crease linearly with density as expected for Debye screening. 



The Ewald convention with 'tinfoil" boundary conditions 
has also been simulated by performing additional integration 
over the mean field E. In this case, all curves collapse onto 
the /f, << L Maxwell result for all temperatures. To reduce 
lattice artifacts at short distances, a dynamic subtraction (that 
will be discussed in a subsequent section) was employed for 
the curves of Fig.|3] 



2. Screening 

One of the characteristic features of systems containing mo- 
bile charges is the phenomenon of screening; the effective in- 
teractions decay exponentially with the inverse Debye screen- 
ing length = ne^/eofesT'^, where n denotes the density 
of the charge carriers. At high enough temperatures the static 
charge-charge structure factor of a screened system has the 
form 



S{q) 



(15) 



In order to test whether our algorithm reproduces correctly 
this essential property of charged systems, we consider a glob- 
ally neutral polyelectrolyte composed of positive and negative 
particles (together with a short ranged Lennard- Jones poten- 
tial) and equilibrate it at the temperature T = 1.25T*. We 
then measure the static structure factor 

S{q) ^ {s{q)s{~q)), s{q) ^ ^=^e^exp{ir^q) (16) 



and compare its form to Eq. il5i in Fig.|4]by plotting 1 / S{q) 
as a function of 1/q^. As can be seen, the long wavelength 
data falls along a straight line, whose slope varies with particle 
density. The solid lines have slope k^L'^/Att'^ and are in good 
agreement with the numerics. 



by unity when measuring density correlations. We also ex- 
amine correlations in the transverse electric field (that part of 
the field for which the discrete divergence, Eq. (|8jl is zero). 
The correlation functions for a single wavevector are shown 
in Fig. 13 a). All modes decay exponentially; time is measured 
in "particle sweeps"; in one time unit we attempt to update 
each particle once together with a variable number of plaque- 
ttes (see below). 

By fitting the correlation functions to exponentials, we ex- 
tract decay rates for selected modes and plot the resulting dis- 
persion relations in Fig. |3b). As can be seen, the density 
modes (□) are diffusive at small wavevectors, lo (x q^, but 
the longitudinal (a) and transverse (0) modes of the electric 
field have a gap at q ~ 0, u! oc ujq + q'^ . This behavior can be 
compared to the continuum equations describing a gas of free 
charge carriers (see also ref. 9). The mean current is driven by 
the electric field and the charge density gradient. 



J = ne^aE - DV p. 



(17) 



n is the number density, D a diffusion coefficient that is re- 
lated to the mobility a via the Einstein relation D — ksTa. 
We insert this relation into the diffusion equation ( I12t to study 
the dispersions of the longitudinal and transverse field compo- 
nents. 

First, we take the curl of Eq. il2\ and Fourier transform to 
find the relation 



(18) 



As observed in Fig. EJb), the transverse modes diffuse with 
a diffusion constant that is determined by a characteristic 
relaxation time scale related to the frequency of updates of 
the transverse field components. In the absence of charges, 
the relaxation time ^ as g ^ 0, but the presence of a 
charge density induces a gap luq ~ Dne^ /ksT — Dk^ in 
the spectrum. The origin of this gap is the screening of the 
interactions on a scale given by the inverse Debye screening 
length K^. 

For the longitudinal field components, we find that 



{iu + Dq^ + ne^a/eo)q • E = 0. 



(19) 
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FIG. 5: (a) Relaxation of particle (thick solid line), charge (thin 
dashed line) and transverse field (thick dashed line) correlations as 
well as exponential fits. System size L = 20 a, p — 0.2 a^'^, 
T — 1.25r*, mode (1,1,1). (b) Dispersion relations for particle 
(□), charge (a) and transverse field (0). Also shown is the density 
dispersion relation for a neutral system (■) with otherwise identical 
parameters. 



Again we find a gap in the spectrum with the same frequency 
ujQ as for the transverse modes. We have also computed nu- 
merically the relaxation of the field correlations (E(t)E(O)) in 
real space, which corresponds to the q = mode in Fig.|5Ib). 
This frequency coincides with ljq. 

Fig. |5] provides the basic demonstration of the efficiency 
of the algorithm. The slowest relaxation time of the system 
corresponds to particle diffusion and scales quadratically with 
system size; any meaningful statistical average must equili- 
brate the sytem over that timescale. The charge and transverse 
fluctuations relax much faster than the density fluctuations, 
and in Monte-Carlo one is free to tune the "diffusion constant" 
(the slope of the straight line) for the transverse field by ad- 
justing the relative frequency of plaquette to particle updates. 
If we update the transverse field very often then we converge 
to the classic algorithms with instantaneous propation of the 
electric field. However, even very rare transverse updates will 
still equilibrate the transverse field components. 

For Fig.|5l we have chosen a high frequency of transverse 
mode updates, where we update typically 2/3 of all plaquettes 
in one particle sweep. Nevertheless, we estimated only about 
15% CPU time is needed for this part of the simulation, all 
other time was consumed by the particle updates (essentially 
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FIG. 6: Dispersion relations in a simulation with system size L = 
20 a, p = 0.01 a"^, T = 0.125r* for particles (□), longitudinal 
(A) and transverse fields {()). For the solid symbols the transverse 
fields were updated less frequently (see text). 



charge interpolation plus the Hamiltonian path). We dropped 
the frequency of plaquette updates by a factor of 10 and ob- 
served no change in the relaxation of the density and charge 
correlations. This result indicates that in a reasonably dense 
system, the particle motion alone is akeady performing a sub- 
stantial part of the integration over the transverse modes. Our 
simple analytic theory predicts exactly this: even in the limit 
of slow propagation of the transverse field the longitudinal dy- 
namics equilibrate. 

As an additional illustration of the efficiency of the algo- 
rithm, we replaced the charged particles of Fig.|5Jb) with neu- 
tral particles and repeated the simulation with otherwise iden- 
tical parameters. The particle dispersion (■) shows that the 
neutral particles diffuse about two times faster, since they are 
not slowed down by the electric field. We conclude that the 
price for simulating this dense charged system is an overhead 
of only about a factor of 2 in particle sweeps (though clearly in 
a neutral system one does not need charge interpolation which 
is relatively costly in cpu time). 

It is clear that the algorithm exhibits strict 0{N) scaling (at 
least away from critical points), since the work for moving a 
particle is purely local and consists essentially in the charge 
interpolation. One might object that the algorithm is surely 
inefficient at very low densities: the effort needed to update 
the plaquettes must become overwhelming. To characterize 
the speed of our algorithm in such an unfavorable situation, 
we contrast the high temperature, high density simulation of 
Fig. |5l with a simulation at a lower temperature T — 0.125T* 
and density p = 0.01 a~'^ in Fig.|S] It is now necessary to 
perform more plaquette updates than before; most particles 
have condensed into pairs that move slowly and cannot easily 
integrate over the transverse field. For the open symbols in 
Fig.|5]we update each plaquette once per particle sweep and 
achieve a very high relaxation rate for the transverse field. The 
computational effort for the transverse updates is about 70% 
of the total runtime in this case. We dropped the frequency of 
transverse updates by a factor of two (filled symbols), but the 
longitudinal and transverse relaxation rates only dropped by 
about 10%. 
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FIG. 7: Dispersion relations for particle (□), charge (a) and trans- 
verse field (0) for simple dipoles. System size L = 20 a, p = 
0.2 a-^ T = 1.25r*, k = SttT*. 




time 

FIG. 8: Squared deviation from Gauss' law/site as a function of time 
in a simulation of size L — 20 a, where the number of particles 
= (solid), iV = 800 (dashed) and = 1600 (long dashed). 



We should note that we are using a very simple Monte- 
Carlo procedure in this low density system, which chooses 
plaquettes randomly in the simulation cell. One can easily 
bias the choice of the plaquette to be in the neighborhood of 
a particle without violating detailed balance. Initial tests on 
our lattice implementation^ have shown substantial increase 
in efficiency at very low densities (volume fraction of 10^^) 
with such a choice. We therefore anticipate that simple, local 
optimizations will greatly improve the performance of the al- 
gorithm for very dilute systems at low T, but leave a detailed 
exploration of this point to future work. 

2. Dynamics of polar molecules 

In the previous paragraph, we argued that relaxation of the 
simple Coulomb gas is aided by screening due to the gap in 
the relaxation spectrum. It is therefore interesting to ask how 
the algorithm performs on a non-screening system. To this 
end, we create polar molecules by coupling two particles with 
opposite charge with a harmonic spring of the form Vb [r) = 
k{r — ro)^. In Fig.^l we show the dispersion relations for this 
polar system. 

Comparing with Fig.|5l we observe that the spectrum for the 
transverse modes (0) no longer has a gap, but rather extrap- 
olates to zero at (7 = 0. However, the dispersion continues to 
be diffusive, and the transverse modes can be relaxed at any 
desired rate by varying the number of plaquette updates per 
sweep. The longitudinal (a) spectrum continues to exhibit a 
gap, but its magnitude is now determined by the spring stiff- 
ness governing the internal vibrations of the molecule. The 
efficiency of the algorithm is independent of the presence of 
screening, and the relaxation of the transverse modes is not 
impeded by the longitudinal field. 

C. Numerical stability 

Another appealing feature of the algorithm besides its sim- 
plicy is its numerical stability. Every upate of the field on a 



link, be it through particle motion or plaquette updates, in- 
volves a numerical error of order the machine precision (typi- 
cally 0(10^^^)). Although this error is cumulative over time, 
it increases so slowly that the accuracy should remain accept- 
able even for the longest run times. To investigate the stability 
of our code, we averaged the deviation of the field configura- 
tion from Gauss's law Eq. (|8|l, 5^ ~ [a^ Ei^j — ei/e^f' 
over the entire simulation cell and plot this quantity as a func- 
tion of time in Fig. |8] Time was again measured in Monte- 
Carlo sweeps. We find a diffusive growth for the error/site, 
jl? = const. X lO^^^t. The constant is smallest for an 
empty system and increases slowly with the number of parti- 
cles, but is always less than one. Since the equilibration time 
is of order 0{l?\ even the largest systems equilibrate with 
acceptable errors. 



IV. DYNAMIC SUBTRACTION OF LATTICE ARTEFACTS 

Many conventional fast electrostatic methods (and the 
present algorithm) interpolate the point charges onto an elec- 
tric grid in order to exploit the fast Fourier transform (FFT)li 
or multigrid algorithms'^. Such an interpolation introduces 
a set of artefacts that have to be controlled: at short scales, 
the interaction between the charge clouds will in general dif- 
fer from that of point charges; furthermore, the self energy 
Usei f of an interpolated charge distribution pint (r) becomes a 
function of its position with respect to the mesh. In general. 

Itself \j d^rd^r' G{V - r')p,nt{Y)p^nt{Y') 

where Smt (q) is the structure factor of the interpolation of the 
charge to the mesh. For the Coulomb interaction, the Green 
functions are GcouK?") ^ l/r and Gcou\{q) ~ l/g^ in real 
and Fourier space, respectively, so that the integral diverges 
linearly for g — > cx). This leads to a periodic potential which 
has an amplitude AUg of order ©(e^/eoa) for mesh size a. 
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Use of a finer mesh leads to worse results with the simplest 
interpolation schemes. The periodic potential tends to "trap" 
the particles in the center of the grid cubes where the self- 
energy is lowest. 

The conventional remedy to the situation consists simply in 
spreading out the charge over a wider range of sites, thereby 
reducing the self-energy artifacts. This is typically done with 
a convolution step that distributes the charges over several 
hundred sites using Gaussians'^ after interpolation to the lat- 
tice. These convolutions are easy to perform in Fourier space, 
but are complicated and time consuming in real space. They 
would also drastically reduce the speed of our algorithm, since 
the Hamiltonian path would have to visit many more sites in 
order to properly update the field configuration. 

Here, we introduce a dynamic correction that leads to an 
effective convolution of the charge distribution with minimal 
computational overhead. For simplicity, we initially work in 
a continuum picture before discretizing and consider the un- 
constrained scalar functional 



pip 



(21) 



Minimizing this functional with respect to the field ip leads to 
the Yukawa equation 



(V^ - H^^JY = -p/to- 



(22) 



Substituting the solution i/jy into Eq. i21\ . we find a negative 
Yukawa Green function. The negative sign of the Green func- 
tion is a feature of the scalar energy functional that we exploit 
here. By coupling the charge density p in our simulation to the 
scalar functional Eq. (12 1> . we generate the partition function 



= / V'ipe' 



and the particles interact with the Yukawa potential 

6162 



Vvir) 



(23) 



(24) 



The essence of our dynamic subtraction algorithm is to cou- 
ple the charges simultaneously to the unconstrained functional 
Eq. (12 1> as well as the constrained vector functional, Eq. ([Q. 
The total partition function then reads 



Z(r) = Zcoulornb{Y) X Zy (t) X COnSt. 



(25) 



and implies an effective interaction between two charges ei 
and 62 given by 



(26) 



At large separations we find the Coulomb interaction; at small 
distances the potential has been regularized. The correspond- 
ing Green's function for ea.( l26> is the sum of the Green's 
functions for Coulomb and Yukawa interaction. 



1 



G(q) = -7 



1 



+ /i^ q^ (q^ + /i^ 



(27) 



H 

CQ 

I 




FIG. 9: Pair potential for combined Coulomb/Yukawa interaction at 
T = 0.0625r*, L ^ 10a, p = 0.75a"\ The top curve (□) 
shows the full combination in agreement with V{r) = const — (1 — 
exTp[— pr]) / ineor ~ r'^ /3eoL^) (solid line). In the lower curve (A), 
the Yukawa interaction was removed for distances up to r < 2.5a. 
For the curves shown in Fig.|3|the Yukawa interaction was removed 
up to L/2. 



By inserting eq. ( I27> into Eq. i20\ . we find that the self-energy 
is now finite for a ^ 0. We can interpret the factor that multi- 
plies the bare Coulomb Green function on the right of ea.( l27t 
as a structure factor Sconv{<i) ~ /^^/('Z^ + u^) describing a 
convolution of the original charge density. In real space, this 
convolution function is 



d^q 

-oo 



.(q) 



l/2g-i:qr 



l^' Kijpr) 
27r2 r '■ 



(28) 



where Ki is a modified Bessel function. 

Use of the scheme Eq. i26i allows us to construct an ef- 
ficient algorithm, since it requires only the introduction of 
a scalar field with Monte-Carlo moves and coupling to the 
charge density according to Eq. ( I2H . The additional overhead 
of the dynamics of this supplementary field is small, since 
we reuse the same values of the interpolated charges for the 
vector and scalar fields. The run time complexity of the dy- 
namic correction is independent of p. Small values of p give 
maximum smoothing, but also strongly alter the interaction at 
small distances. In general, one therefore wants to remove the 
Yukawa interaction on a scale of 1/p by adding —Vyir) to 
the interaction. This can be conveniently done together with 
all the other short range (e. g. van der Waals) interactions of 
the system. The choice of p therefore gives us a freely tunable 
tradeoff between precision and efficiency. This algorithm re- 
minds us of early attempts at creating finite field theories in 
the 1940's^§ using different particles to cancel divergences in 
physical quantities. Here it is the ground state energies of two 
Gaussian field theories (but with scalar and vector symme- 
tries) which have compensating divergences in the continuum 
limit, leading to a finite theory when they are combined. 

We have implemented the above method in our off-lattice 
Monte-Carlo simulation to test its validity and efficiency. To 
this end, Monte-Carlo moves for the second auxiliary field ip 
were introduced to generate the effective interaction Eq. ( I26> . 
We first verified this interaction by extracting the pair po- 
tential from the distribution of separation P{r) of a pair of 
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FIG. 10: Self energy along the diagonal of a grid cube from comer 
(0, 0, 0) to center (0.5, 0.5, 0.5) for the pure Coulomb interaction 
(lowest curve) and the combined Coulomb/Yukawa interaction for 
H = 3a~^, fj, = la~^, and /i = 0.5a~^ (highest curve). Charges 
were interpolated using B-splines, and the energies were calculated 
using the lattice Green function Eq. <29> . 



charges as in Sec. 1111 A 1 1 The top curve in Fig. |5] shows 
the combined Yukawa and Coulomb interaction for a value 
of n = 1 in the low temperature limit of Maxwell boundary 
conditions. The Yukawa field was updated much more rarely 
than the transverse electric field since the massive field relaxes 
rapidly. The curve agrees well with the analytic expression 
V{r) — const — (1 — exp[— /xr])/47reor — r'^ /3eoL^) ex- 
cept at small distances where the interpolated charges overlap. 
The Yukawa potential can be subtracted Eq. i24l explicitly at 
short distances. The lower curve shows a case where this was 
done for r < 2.5a; the numerical points agree well with the 
Coulomb potential. Note that this correction at small scales 
should in fact not be neccessary in realistic applications that 
use a grid size a that is much smaller than the LJ particle radii. 

We now have to show that the imposition of an additional 
Yukawa interaction also reduces the self-energy artifacts. The 
most dramatic manifestation of these artifacts occurs through 
a trapping of particles in the cell center due to the (periodic) 
energy baiTier AUs = J7cornor - C^ccntcr, where C/contor and 
C^cornei are the self energies of a particle located at the center 
or corner of a grid cube. The self energy of a particle whose 
charge is interpolated onto the cubic grid with the B-spline 
interpolation Eq. can be calculated using the lattice Green 
function for the Yukawa interaction. 



(27r)3 6 — 2 cos g^; — 2 cos Qy ~ 2 cos qz + jJ. 



We evaluated Eq. ( I29> numerically using Mathematica. 
Fig. ^1 shows the self-energy along the diagonal of a grid 
cube for for the combined Coulomb/Yukawa interaction for 
several values of /i. The amplitude of this potential is largest 
for the pure Coulomb case and is drastically reduced when 
fi = 0.5a-i 

In our Monte-Carlo simulations, AUg has been extracted 
from the ratio of probabilities to find the particle in the cen- 
ter or corner, respectively. In Fig. ^2 we show the effective 




FIG. 1 1 : Reduction of the self energy barrier leading to trapping of 
particles in cell centers. Numerical points correspond to tempera- 
tures T = 0.0625r* (□) and T = 0.125r* (A). The probabili- 
ties were measured by counting particles in cubic regions centered 
in center and comer of size 0.4'^a'^ = 0.0640"^. The lines show the 
expected reduction of the energy barrier based on a semi-analytic 
evaluation of the self-energy using the scheme Eq. <26> (solid line) 
and Eq. <31> (dashed line). 



self-energy barriers obtained in this way as a function of /i. 
As expected from Fig. ^| the barrier drops rapidly with de- 
creasing fi, and almost complete elimination is achieved at a 
value of /i = 0.5a^^. The solid line in Fig.^Jwas obtained 
by computing the self energy difference of particles located 
in the center and corner of the grid cube as shown in Fig. 1101 
The numerical points agree well with this result. 

From the numerical integration of Eq. ( I29l l, we find that the 
amplitude of the periodic energy barrier varies as AUs{n) = 
2 X 10^^ (^a)^T* for small fi. Even better smoothing should 
be possible if one uses a mixture of Green functions that de- 
cays even faster for large wavevectors than Eq. ( I27l l: Let us 
consider 



G(q) 



2m^ 



q2(q2 -I- //2)(q2 _^ 2^^) ' 
The corresponding real space potential is 



^ 6162 

= 1 (1 



(30) 



(31) 



which can be generated by coupling the charges to three 
fields, respectively a constrained vector, a scalar and an un- 
constrained vector. This unconstrained vector functional reads 



eoE(r)2 (eodivE(r)-p(r))2 



2/^2 



(32) 



which leads to a Yukawa Green function with a positive pref- 
actor. The dashed line in Fig. shows the energy barrier 
as a function of /i for this scheme. The inhomogeneity in 
the self energy is smaller than when using only the scalar 
Yukawa functional, and the barrier now varies as AL/s(/i) = 
1.5 X 10~^ {iia)'^T* for small /i. Although we have not im- 
plemented this full scheme in our Monte-Carlo simulation, it 
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appears to be the strategy of choice when optimizing the algo- 
rithm for truly atomistic, high accuracy simulations. We note 
that the sucess of the method requires that both the vector and 
scalar functionals implement equivalent discretizations of the 
Laplacian operator. In our code this is the simplest 7-point 
difference scheme. 

We would like to emphasize that the present strategy is gen- 
eral and can also be used in combination with other Coulomb 
algorithms. Methods relying on solutions of the Poisson equa- 
tion could simulate an additional Yukawa equation and add 
the resulting potentials or forces. This seems to be rather 
attractive, since the charge spreading step amounts to a sig- 
nificant portion of the total computational effort in multigrid 
formulations^. 



V. CONCLUSIONS 

We have implemented a Monte-Carlo algorithm for charged 
particles moving in the continuum. Instead of minimizing the 
electrostatic energy, the algorithm uses particle moves in com- 
bination with constrained auxiliary field updates whose dy- 
namics are inspired by the full Maxwell equations. Unlike 
in electrodynamics, however, the field propagates diffusively. 
We showed that the algorithm efficiently equilibrates charged 
systems such as electrolytes and polar fluids. It is most effi- 
cient at higher densities and temperatures, where the particle 
motion alone can perform a substantial part of the integration 
over the transverse modes of the electric field. Even at lower 
temperatures, the cost for the transverse integration remained 
acceptable. 

We also introduced an approach to treat lattice artifacts that 
arise from the charge interpolation onto a grid. By introducing 
an additional negative Yukawa interaction, the periodic poten- 
tial arising from the self-energy of the particles can be reduced 
to any desired accuracy. The scheme is completely general, 
and other Coulomb algorithms could also benefit from the 
method. 

The principal advantages of the present approach are 

• True 0{N) scaling with the number of particles in the 
system. No such algorithm is available for Monte- 
Carlo simulations at present. Many problems involving 
electrostatics that could previously only be treated ef- 
ficiently with molecular dynamics can now be treated 
with Monte-Carlo, which is often easier to program and 
faster due to larger time steps. 

• Purely local computations are required to advance parti- 
cle and field degrees of freedom, making the algorithm 
an ideal candidate for parallelization on large compute 
clusters or supercomputers. 

In the present form, the algorithm can treat coarse-grained 
models of e. g. polyelectrolytes or membranes. It now has 
to pass more stringent tests with realistic model systems and 
truly atomistic simulations. 
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APPENDIX A: BOUNDARY CONDITIONS FOR CHARGED 
PERIODIC SYSTEMS 

We first summarize the boundary conditions employed 
for Ewald summation and then discuss their relation to the 
Maxwell boundary conditions, which are the natural choice 
for our algorithm. The Ewald formula for the electrostatic en- 
ergy of periodic a system is 



U ^ eiCj^lperi'Ti - Vj) + Uself + Udip (Al) 



The self-energy, Useif — '^/v^ X^i , is constant and 
therefore irrelevant for Monte-Carlo dynamics. The periodic 
Ewald potential ipper consists of sums in real and Fourier 
space: 



47r ^ exp(-G'V4K2 + iG • r) 



TT 



E 

R 



Q2 

erfc(K|r + R| 



Rl 



(A2) 



The result of the sum is independent of the range parameter, k. 
At short distances, this expression can be expanded'^, giving 



■0 = const - 



1 



27rr^ 



(A3) 



The second term is of the same order in 1/L as the dipole 
energy 



1 



2eo(2e' + 1)L3 



E 



(A4) 



This energy depends on the dielectric constant e' of the sur- 
rounding medium. Two common choices for the dielectric 
constant are e' = 1 ("vacuum") and e' = oo ("tinfoil"). 

A term comparable to Udip arises when using Maxwell 
boundary conditions: In SecHlJ we saw that the field dynam- 
ics obey a Maxwell equation. 



i9E 

— = -J/eo + curlB//io- 
ot 



(A5) 



Integrating this equation over space and time, we find the re- 
lation between the q = component of the electric field and 
the dipole moment, d 



(A6) 
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since the integral of the curl of a periodic function is zero. The 
contribution of E to the electrostatic energy is 



fo 
2 



(A7) 



The longitudinal components of E for q 7^ are given by the 
gradient of iljper/4:TTeo. 

Note that there is an important difference between the con- 
ventional dipole moment of the simulation cell, ^ e^ri and 
d: A charge can wind about the periodic system giving equiv- 
alent configurations with different energies. These two def- 
initions of the dipole moment are only the same modulo a 
Bravais lattice vectoriiS: 



E 



|e|L{Z, m, n}. 



(A8) 



for arbitrary integers /, m, n. 

Consider now a pair of particles with opposite charges sep- 
arated by r, simulated by our algorithm. If we sum over all 
equivalent configurations we find the relative statistical weight 
of a given state 



-e^{r+L{Lm,n})^/2kBTtoL^ 



(A9) 



At low temperatures and with r small the sum Eq. (IA9> is 
dominated by a single minimum energy term with I = rn = 
n — 0, thus 

Zo = e-'=''-V2fcBT<:oL=' foj. fc^j^ ^ eV4eoi (AlO) 



The corresponding free energy is 

F = -fcTlnZo = e^r^/2eoL^. 



(All) 



We now combine the low temperature result with the short 
distance expansion of the Ewald potential for two particles by 
adding Eq. (lAlU to Eq. (IA3> . We find 



V{r) 



f 1 



V47reor 6eoL^ J 2eoL^ 



(A12) 



We thus expect a net quadratic correction e^r^ /3eQL^ to the 
l/r interaction as in Eq. ( I14> . 

In the opposite limit of high temperatures, the charges are 
mobile and the current can wind around the simulation cell 
many times. In this limit 



Zq const, for ksT e^/4eoL. 



(A13) 



There are no dipole contributions in this limit, and the poten- 
tial is that of Eq. ( fA3t . 

Note that the above arguments hold for free charges that 
allow the current to flow freely; for bound charges or dipoles 
the current does not lead to multiple inequivalent states for the 
field E. It is possible to impose the tinfoil boundary condition 
in our algorithm by introducing a third Monte-Carlo move that 
integrates over the total electric field E (see also Ref.Q3). 
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